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ABSTRACT 
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per  km.  Several  new  parameters  for  use  in  mode  search  are  deduced  from  the  results  and 
some  old  ones  are  verified.  Future  studies  in  waveguide  mode  propagation  theory 
pertaining  to  atmospheric  ducts  may  benefit  from  this  work.  An  improved  mode  search 
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I.  INTRODUCTION 

The  M-Layer  program  developed  by  NOSC  (Naval  Ocean  System  Center)  was 
documented  by  Yeoh  [Ref.  1]  at  NPS.  It  was  later  extensively  revised  by  Lee  and  Han 
[Ref.  2]  for  improved  accuracy  and  speed.  The  new  FORTRAN  code  is  now  identified 
as  the  NPS  version  [Ref.  3]  under  the  auspices  of  NRaD  (Naval  Command,  Control  and 
Ocean  Surveillance  Center,  RDT&E  Division),  the  current  name  for  NOSC.  In  this 
version,  the  logical  structures  and  numerical  algorithms  for  following  the  constant  phase 
lines  of  the  mode  function  and  for  computing  the  mode  locations  were  almost  completely 
rewritten.  The  overall  plan  to  first  partition  the  search  region  into  rectangular  areas  and 
then  circulate  around  these  rectangles  to  search  for  the  modes,  i.e.,  the  mode  search 
protocol,  was  left  unchanged.  Lee  and  Han  [Ref.  2]  recommended  that  a  more  direct 
mode  search  protocol  should  be  devised  to  locate  the  modes  more  expediently.  If  possible, 
such  an  approach  should  locate  the  modes  in  the  order  of  ascending  range  attenuation 
rates. 

To  design  a  more  efficient  mode  search  protocol,  a  better  understanding  of  the  mode 
function  beyond  the  known  locations  of  the  modes  is  required.  In  this  thesis,  the 
analytical  property  of  the  mode  function  is  investigated.  Programs  are  written  to  track  and 
plot  the  constant  phase  lines  along  which  the  mode  function  is  either  real  or  purely 
imaginary.  From  the  results,  a  new  strategy  to  locate  the  modes  is  recommended. 


In  the  remainder  of  this  chapter,  the  background  theory  and  some  of  the  notations 
used  in  this  thesis  are  introduced.  The  basic  theory  presented  in  Sections  A  and  B  follows 
Ref.  3  closely.  The  results  of  this  research  are  presented  in  Chapter  II.  In  Chapter  IE,  a 
new  mode  search  protocol  is  proposed  based  on  the  findings  of  this  thesis.  The  relevant 
parameters  for  use  in  this  new  approach  are  also  discussed. 

A.      THE  WAVEGUIDE  MODE  THEORY  OF  PROPAGATION 

Trapping  of  electromagnetic  (EM)  waves  in  the  modes  supported  by  a  duct  is  the 
dominating  factor  in  over-the-horizon  propagation.  The  computer  program  M-Layer 
searches  for  these  modes  and  computes  the  electric  field  from  the  Hertz  vector.  Assume 
a  vertical  electric  dipole 

J  =  4nzt>(x)&(y)6(z-zT) 

at  a  height  z=Zy  above  the  surface  of  the  earth.  Following  Freehafer  [Ref.  4],  the  Hertz 
vector  of  the  EM  fields  can  be  written,  in  the  cylindrical  coordinates  (p,  <|>,  z),  as  a  sum 
of  contributions  from  individual  modes  [Ref.  1]: 

H(p,z)  =  -*;£  ff0P)(PmP)Sm(Zr)Sjz)  (1) 

m 

where  Pm  is  the  wavenumber  of  the  m-th  mode  and  is  independent  of  the  coordinates; 
Zj  is  the  height  of  the  transmitter  above  the  ground,  Hq'  >  is  the  Hankel  function  of  the 
second  kind,  which  represents  an  outgoing  wave  in  the  radial  direction  when  the  time 


dependence  eJ"1  is  adopted,  and  gm  is  the  height-gain  function  of  the  m-th  mode, 
normalized  so  that  the  integral  of  its  square  over  all  heights  equals  unity  when  either  an 
electric  or  a  magnetic  dipole  source  is  used.  The  height-gain  function  satisfies  the 
equation 


\dz 


gj®=0  (2) 


where  k  is  the  free- space  wavenumber,  m(z)  is  the  modified  index  of  refraction,  and 
m  (z)  is  approximated  with  a  continuous,  piecewise  linear  profile  having  I  layers: 


m2(z)=m,2 +0,(2-2,)  z^zzz^,     lzizl.  O) 

In  practice,  m(z)  deviates  only  slightly  from  unity  in  the  troposphere.  The  modified 
refractivity  M(z)  =  [m(z)-l]xlO  ,  as  shown  in  Fig.  1,  is  commonly  used.  The  slope  of 
M(z)  is  approximately  cc-xlO  /2. 

In  the  i-th  layer,  the  height-gain  function  is  given  by 

where  q^  is  a  dimensionless  linear  function  of  height  z  with  k,  P  m-,  and  or-  as 
parameters: 
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Figure  1.  Typical  modified  refractivity  M(z)  profile. 
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The  functions  kjCq^)  and  k2(qmj)  are  given  in  terms  of  the  Airy  function  Ai: 

kl(qJ=-j2(l2)ll6Ai(-qmi^2^)  (6) 

and 

^(0=2(12)1/6Ai(-0.  (7) 

For  z  >  zj+j,  irr(z)  is  extended  continuously  upward  at  a  constant  slope 
commensurate  with  the  effective  radius  of  the  earth.  This  slope  equals  2.36x10  (m~  ) 
for  the  four-thirds  effective  earth  radius  model.  Thus,  the  height-gain  function  gm  is  again 
given  by  Eqs.  (4)  through  (7)  for  z  >  Zj+j,  with  the  constant  Cj+j  set  to  e^  n'  to 
represent  an  upward  going  wave  as  z  becomes  large.  Beneath  the  "flattened"  earth  surface 
where  z  <  Zj  =0,  m  (z)  is  deduced  from  the  dielectric  constant  and  the  conductivity  of 
sea  water,  which  are  set  to  80.8869  and  4.64  (Si/m)  respectively.  Here  the  Hertz  vector 
is  represented  by  a  downward  propagating  plane  wave  in  this  region. 

The  conditions  that  the  tangential  components  of  the  electric  and  magnetic  fields 
are  continuous  across  the  layer  boundaries  determine  the  coefficients  C-.  Starting  from 
the  top  layer  down,  or  "integration-down"  [Ref.  1],  every  Cj  is  uniquely  determined  by 
Cj+^  at  z  =  Zj+1.  Thus,  a  set  of  all  C-  coefficients  is  determined  without  considering  the 
boundary  conditions  at  z  =  Zj.  On  the  other  hand,  starting  from  the  lowest  boundary  at 


z  =  Zj  up  to  z  =  Zj,  or  "integration-up",  a  second  set  of  C^  coefficients  is  determined 
without  applying  the  boundary  conditions  at  z  =  Zj+j.  That  these  two  sets  of  Cj 
coefficients  are  identical  is  a  manifestation  that  Pm  is  the  wavenumber  of  a  mode. 
Mathematically,  this  consistency  condition,  sometimes  called  the  guidance  condition 
[Ref.  5],  can  be  expressed  as  the  vanishing  of  a  determinant  D  called  the  mode  function 
[Ref.  6].  The  elements  of  this  determinant  consist  of  linear  combinations  of  kjCq^)  and 
k2(qrm)  and  their  derivatives  with  respect  to  height  at  the  boundaries.  The  M-Layer 
program  searches  for  the  zeros  of  the  mode  function  to  find  the  wavenumbers.  Once  a 
wavenumber  is  obtained,  the  height-gain  function  of  the  corresponding  mode  can  be 
computed.  Note  that  the  normalization  condition  on  gm  and  the  boundary  conditions, 
together  with  the  C^  coefficients,  determine  the  B-  coefficients  to  within  an  overall  sign. 
This  sign  need  not  be  resolved  because  only  products  in  the  form  of  gm(z)gm(zj)  from 
each  mode  contribute  to  the  Hertz  vector. 

Since  the  mode  function  D  depends  on  the  wavenumber  Pm  only  through  the  values 
of  q^  at  the  layer  boundaries,  it  is  more  convenient  to  consider  D  as  a  function  of  the 
variable  q  given  by 
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a, 


•'-£, 


(8) 


and  to  search  for  the  zeros  of  D(q)  in  the  complex  q  plane.  The  m-th  zero  is  designated 
as  qm  and  is  called  a  q-eigenvalue  and  Pm  is  then  deduced  from  qm  by  inverting  Eq.  (8) 


for  p.  In  the  NPS  version,  q     is  ordered  in  ascending  attenuation  rate  of  the  mode,  which 
is  approximately  proportional  to  the  imaginary  part  of  Pm. 

B.      MODE  SEARCH  PRINCIPLE 

M-Layer  searches  for  modes  which  have  attenuation  rates  below  a  specified  value. 
At  ground  ranges  more  than  several  wavelengths  away,  this  attenuation  rate  is 
proportional  to  the  imaginary  part  of  the  wavenumber  pm  as  can  be  seen  from  Eq.  (1). 
In  the  upper  complex  q  plane,  for  values  of  q  such  that 


it2 


«1. 


(9) 


mj-P/k  is  approximately  proportional  to  q  because  the  absolute  value  of  m^  and,  hence, 
that  of  p/k,  are  both  nearly  unity.  Thus,  the  limit  on  attenuation  rate  imposed  on  the 
imaginary  part  of  p  can  be  translated  approximately  into  an  upper  bound  on  the 
imaginary  part  of  q  which  defines  a  strip  in  the  complex  q  plane  called  the  search  region. 
This  region  is  covered  with  layers  of  identical  square  meshes  whose  sides  are  parallel  to 
the  imaginary  q  axis  and  have  a  length  equivalent  to  an  attenuation  rate  of  1/32  dB    per 


Since  the  absorption  is  small  in  air,  the  modified  index  of  refraction  m(z),  and 
hence  8|,  are  considered  as  real  quantities  in  the  following  discussions.  In  the  actual 
FORTRAN  code,  they  are  declared  as  complex  variables. 

This  is  the  default  value  of  the  NPS  version  which  can  be  adjusted  through  editing 
the  variable  "dmesh"  in  an  ASCII  input  file. 


3 

kilometer.  The  lower  edges  of  the  meshes  in  the  lowest  layer  sit  along  the  real  q  axis  . 
The  meshes  in  the  layer  just  below  the  top  one  contain  the  upper  bound  of  the  search 
region,  with  the  top  layer  providing  some  allowance  for  numerical  inaccuracy.  The 
program  tests  each  mesh  square  for  the  presence  of  zeros  of  the  mode  function  D(q). 

To  search  through  all  the  meshes,  the  program  first  divides  the  mesh  covering  the 
search  region  into  "contour  rectangles"  with  equally  spaced  vertical  lines  parallel  to  the 
imaginary  q  axis.  These  vertical  lines  contain  the  edges  of  stacked  square  meshes  and  are 
separated  by  a  distance  160  times  the  side  of  a  mesh.  The  search  commences  at  the  top 
left  corner  and  moves  counterclockwise  around  each  "contour  rectangle"  and  begins  with 
the  one  whose  left  edge  is  defined  by 


Re 


(1X 


-2/3 


=/&  mj-m^  -2xlO-6l2e(M1-Afinin)=2xlO-6dl 


(10) 


where  mn^n  is  the  minimum  value  of  the  modified  index  of  refraction  profile  and  dmmm 
is  shown  in  Fig.  1.  The  justification  for  this  choice  to  begin  the  mode  search  is  as 
follows:  From  Eq.  (5)  and  optics,  for  a  wave  to  propagate  freely  in  space,  the  dispersion 
relation  of  the  Maxwell  equations  requires  that  k^m^(z)  >  P  m,  assuming  that  both 
quantities  are  real.  Thus,  the  smallest  p  m  to  support  a  trapped  wave  will  occur  when 


3 
The  NOSC  version  extends  the  lower  edge  slightly  below  the  real  q-axis.  This 

causes  problems  in  some  situations  [Ref.  2]. 
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0  just  exceeds  the  minimum  of  m  (z)  .  It  is  not  surprising  that  an  argument  based 
on  optics  works  within  a  waveguide  mode  formulation  which  is  low  frequency  in 
principle.  Lee  [Ref.  7]  has  demonstrated  that  the  earth-flattening  approximation  actually 
links  Mie's  low  frequency  oriented  spherical  harmonics  to  Fock's  high  frequency 
diffraction  theory  through  the  uniform  asymptotics.  The  mere  introduction  of  an  unlimited 
ground  range  into  the  Maxwell  differential  equations  incorporates  the  global  feature  of 
the  radius  of  the  earth  into  the  local  equations. 

After  the  search  over  the  initial  rectangle  is  completed,  the  program  goes  on  to 
search  the  neighboring  rectangle  to  the  left.  If  a  specified  number5  of  consecutive 
rectangles  of  decreasing  real  q  values  have  been  searched  without  turning  up  any  zero  of 
D(q),  the  program  changes  direction  and  starts  to  search  the  rectangles  to  the  right  of  the 
initial  "contour  rectangle"  one  by  one,  with  increasing  real  part  of  q.  After  failing 
consecutively  to  locate  any  q-eigenvalue  again  over  the  specified  number  of  rectangles, 
the  program  assumes  that  no  more  zeros  of  D(q)  exists  in  the  search  region.  The  mode 
search  is  considered  complete  and  the  procedure  is  terminated.  If  the  array  for  storing  the 
q-eigenvalues  is  filled  up  before  the  search  is  completed,  the  search  is  terminated  with 
an  error  message. 


Note  that  assuming  the  same  real  part,  an  imaginary  component  of  fim  reduces  the 
real  part  of  P  m  which  may  enable  the  wave  to  propagate  in  the  z-direction. 

5  This  number  of  consecutive  rectangles  is  an  adjustable  input  variable  "nstop"  in  the 
NPS  version.  The  default  value  is  2. 


The  search  for  zeros  of  D(q)  makes  use  of  the  fact  that  a  real  valued  function 
changes  sign  when  it  crosses  a  simple  zero.  Since  a  zero  of  the  complex  valued  function 
D(q)  is  where  both  its  real  part  and  imaginary  part  vanish,  a  necessary  condition  for  a 
point  qm  to  be  a  zero  is  that  it  is  the  intersection  of  two  curves  defined  by  Im{D(q)}  =  0 
and  Re{D(q)}  =0.  M-Layer  moves  along  each  side  of  a  "contour  rectangle"  while 
searching  for  a  sign  change  in  Im{D(q)}  across  an  edge  of  a  mesh  bordering  the  side  of 
this  rectangle  to  determine  that  a  line  of  Im{D(q)}  =  0  has  been  encountered.  The  search 
then  follows  this  line  into  the  meshes  within  the  "contour  rectangle",  checking  each  mesh 
to  see  if  a  curve  Re{D(q)}  =  0  enters  the  mesh  being  inspected.  All  these  steps  make  use 
of  the  assumption  that  the  zeros  of  D(q)  are  simple.  Once  both  the  curve  Im{D(q)}  =  0 
and  the  curve  Re{D(q)}  =  0  are  found  to  be  present  within  a  mesh,  the  locations  of  their 
possible  intersections  are  estimated. 

The  rule  for  sign  change  becomes  inconclusive  if  a  zero  of  Im{D(q)}  =  0  or  if 
Re{D(q)}  =  0  happens  to  fall  on  a  corner  of  a  mesh,  and  a  remedy  is  required.  In  the 
NPS  version,  whenever  the  real  part  or  the  imaginary  part  of  D(q)  vanishes  on  a  comer 
of  a  mesh,  the  phase  angle  of  D(q)  is  rotated  by  2  radians.  This  maneuver  effectively 
shifts  q  by  a  small  amount  to  resolve  the  sign  ambiguity. 

To  estimate  the  locations  of  zeros  of  D(q)  within  a  mesh,  D(q)  is  assumed  to  be 
well  approximated  in  this  mesh  by  its  four-term  Taylor  series  expansion.  The  unshifted 
values  of  D(q)  at  the  mesh  corners  determine  this  cubic  polynomial  uniquely.  Cardan's 
formulas  [Ref.  8]  are  used  to  locate  the  zeros  of  this  polynomial  before  retaining  only 
those  lying  within  the  mesh. 
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IL  MODE  FUNCTION  IN  THE  COMPLEX  qjj/tj  PLANE 

A.      MODE  LOCATIONS 

M-Layer  searches  for  the  zeros  of  the  mode  function  D(q)  in  the  complex  q  plane 
by  tracking  the  constant  phase  lines  of  D(q)  along  which  the  mode  function  is  either  real 
or  purely  imaginary.  The  zeros  found  are  called  the  q-eigenvalues  as  defined  by  Eq.  (8) 
and  are  denoted  as  qm.  These  eigenvalues  are  saved  in  an  ASCII  file  and  are  utilized  later 
for  height-gain  function  computations.  In  order  to  design  a  strategy  for  locating  these 
zeros,  the  mode  locations  are  plotted. 

It  is  clear  from  Eq.  (8)  that  the  variable  q  varies  with  a^,  the  slope  of  m  (z)  in  the 
lowest  (first)  layer.  Since  Oj  depends  on  how  the  continuous,  piecewise  linear 
approximation  to  m  (z)  is  made,  while  both  m^  and  pm  are,  in  principle,  dependent  only 
on  the  actual  profile,  q  will  vary  strongly  with  a  i  and  fail  to  provide  information  on  the 
wavenumber  of  the  mode  directly.  A  more  suitable  variable  to  use  is,  from  Eq.  (8): 


f4%'-£.  (id 


Since  q  is  given  the  variable  name  of  q^  and  (k/a  p2'3  is  given  the  variable  name  of  tj 
in  the  M-Layer  FORTRAN  code,  this  variable  is  identified  as  q^j/tj  throughout  this 
thesis.  Since  m1  is  real  and  both  it  and  0/k  deviate  from  unity  only  slighdy  for  all  cases 
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investigated,  q^j/tj  is  approximately  2(rrij  -  /fyk).  Thus,  its  imaginary  part  is  directly 
proportional  to  -p,  the  attenuation  rate  of  the  mode.  Furthermore,  qj  j/tj  does  not  depend 
explicitly  on  a^  Removing  the  tj  factor  from  q^j  makes  it  possible  to  compare  results 
from  different  evaporation  duct  profiles  meaningfully. 

Plots  of  the  q-eigenvalues  as  individual  points  in  the  complex  qij/tj  plane  are 
included  in  Appendix  A.  For  the  20  meter  duct,  which  will  be  used  as  the  representative 
case  for  discussions  in  this  thesis,  a  line  linking  all  the  eigenvalues  in  the  order  of 
increasing  attenuation  rate  is  shown  in  Fig.  2.  Because  of  the  long  excursion  between 


Figure  2.  q-eigenvalues  in  the  qjjAj  plane,  connected  in  the  order  of  ascending 
attenuation  rate  (20  m  duct). 
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many  pairs  of  consecutive  modes  in  the  upper  part  of  this  figure,  it  appears  that  the 
intuitive  approach  of  starting  with  the  mode  of  the  lowest  attenuation  and  continually 
looking  for  the  mode  with  the  next  higher  attenuation  rate  may  not  be  the  most  efficient 
strategy.  Furthermore,  there  seems  to  be  a  fork  near  the  middle  of  the  figure  which 
suggests  that  there  can  be  more  than  one  constant  phase  line  passing  through  all  the  zeros 
of  the  mode  function.  These  problems  suggest  that,  in  order  to  design  an  efficient  mode 
search  strategy,  the  behavior  of  the  mode  function  in  the  complex  q^/tj  plane  has  to  be 
investigated  more  thoroughly. 

B.      PHASE  LINE  TRACKING 

The  M-Layer  program  follows  a  line  along  which  Im{D(q)}  =  0  until  a  zero  is 
encountered  at  its  intersection  with  another  line  along  which  Re{D(q)}  =  0.  Since  there 
is  no  way  to  determine  the  phase  of  the  constant  phase  line  linking  two  adjacent  zeros  a 
priori,  following  these  two  types  of  somewhat  arbitrary  but  easy  to  compute  phase  lines 
which  have  constant  phases  of  0  or  ji,  and  tt/2  or  3tt/2,  respectively,  are  still  the  best  way 
to  program  a  computer  to  locate  the  zeros  systematically.  It  is  thus  necessary  to  find  out 
the  actual  distribution  of  these  constant  phase  lines  in  the  complex  q^/tj  plane. 

1.       Downward  Tracking 

In  the  original  design,  along  the  real  q  axis,  M-Layer  divides  the  search  region 
into  "contour  rectangles,"  each  of  which  spans  160  meshes  horizontally.  From  the 
observed  locations  of  the  modes,  it  is  found  desirable  to  track  and  plot  the  constant  phase 
lines  along  Im{D(q)}  =0  and  Re{D(q)}  =0  over  a  little  more  than  six  "contour 
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rectangles",  roughly  three  to  the  right  and  three  to  the  left  of  the  starting  real  q  coordinate 
given  by  Eq.  (10),  which  is  a  variable  called  qtest  in  M-Layer.  The  subroutine  FNDMOD 
and  FZEROX  are  modified  to  track  these  constant  phase  lines.  A  listing  of  these  modified 
routines  for  tracking  the  lines  Im{D(q)}  =  0  beginning  from  the  top  edge  and  moving 
downwards  into  the  search  region  is  included  in  Appendix  B.  Their  flow  charts  are  given 
in  Appendix  C.  Only  minor  changes  are  required  to  track  the  type  of  lines  Re{D(q)}  =  0, 
or  to  track  these  lines  from  the  bottom  of  the  search  region  upwards.  The  program  first 
searches  for  a  sign  change  in  Im{D(q)}  or  in  Re{D(q)}  along  the  top  edge  of  the  search 
region,  starting  at  the  coordinate  -512  mesh  sizes  from  qtest,  until  a  distance  of  1024 
mesh  sizes  is  covered.  When  a  sign  change  is  observed,  a  constant  phase  line  is 
recognized  as  passing  through  the  particular  mesh  square  and  the  program  writes  the 
complex  q  value  of  the  lower-left  comer  of  this  mesh  square  into  an  ASCII  file  identified 
by  whether  the  mode  function  is  real  or  imaginary  along  the  line,  and  the  order  this  line 
is  found.  The  program  then  moves  from  the  top  edge  downwards  to  follow  this  constant 
phase  line  until  it  exits  the  search  region.  The  q  values  of  the  lower-left  comer  of  the 
mesh  squares  along  this  phase  line  are  also  written  into  the  same  ASCII  file  to  be  plotted 
later  using  MATLAB/386.  The  constant  phase  lines  for  the  mode  functions  of  the  2,  4, 
6,  8,  10,  20,  30  and  40  meter  evaporation  ducts  are  obtained  and  plotted.  They  are 
included  in  Appendix  D.  Tracking  and  plotting  these  constant  phase  lines  are  extremely 
time  consuming.  The  CPU  time  required  to  track  the  desired  constant  phase  lines  for  each 
duct  using  an  Intel  80486  based  PC  running  at  a  33  MHz  clock  rate  is  listed  in  Table  1. 
The  time  required  to  plot  those  lines  for  each  duct  using  an  Intel  80386  based  PC  running 
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at  a  16  MHz  clock  rate  is  also  listed. 


TABLE  1 
CPU  Time  for  Tracking  and  Plotting  Constant  Phase  Lines 


duct 
height 

(m) 

tracking 
Im{D(q)}  =  0 

(hr:min:sec) 

tracking 
Re{D(q)}  =  0 

(hr:min:sec) 

plotting 
both  lines 

(hr:min:sec) 

2 

1:47:39 

1:46:23 

0:21:28 

4 

2:47:41 

2:48:50 

0:23:17 

6 

3:43:54 

3:45:11 

0:25:42 

8 

4:25:49 

4:23:41 

0:29:34 

10 

6:08:14 

6:11:03 

0:32:25 

20 

8:26:58 

8:24:22 

0:34:21 

30 

9:08:14 

9:06:45 

0:36:49 

40 

8:39:04 

8:37:54 

0:36:53 

Total 

45:07:33 

45:04:09 

4:01:29 

The  constant  phase  line  plot  for  the  mode  function  of  the  20  m  evaporation 
duct  is  given  in  Fig.  3.  The  solid  lines  represent  those  having  Im{D(q)}  =  0;  the  dotted 
lines  represent  those  having  Re{D(q)}  =  0.  Every  intersection  of  these  two  types  of  phase 
lines  is  the  location  of  a  q-eigenvalue  scaled  by  tjf  thus  indicating  the  presence  of  a 
mode.  Note  that  every  solid  line  intersects  with  a  dotted  line  except  for  many  of  those 
running  from  the  top  through  the  bottom  of  the  search  region.  None  of  the  solid  lines 
intersect  with  other  solid  lines,  neither  do  the  dotted  lines.  This  is  seen  clearly  in  Fig.  4, 
which  magnifies  the  lower  center  portion  of  Fig.  3.  One  feature  that  can  be  recognized 
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Figure  3.  Constant  phase  lines  in  the  qj  j/tj  plane  (20  m  duct). 
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Figure  4.  Magnified  lower  center  portion  of  Fig.  3  showing  no  intersection  of  phase  lines 
of  the  same  type. 
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immediately  in  all  the  plots  in  Appendix  D  is  that  the  starting  real  q  coordinate  given  by 
Eq.  (10),  which  is  marked  by  an  "x"  along  the  top  edge  of  each  plot,  correlates  very 
closely  with  the  mode  of  minimum  attenuation. 

2.       Upward  Tracking 

One  mode  of  low  attenuation  rate  is  missing  from  Fig.  2  compared  to  those 
found  in  Ref.  2.  Several  of  them  are  also  missing  from  those  cases  of  greater  duct  heights 
plotted  in  Appendix  D.  It  is  evident  that  tracking  the  constant  phase  lines  from  the  real 
qj|/tj  axis  upwards  is  necessary.  For  the  10  meter  and  the  20  meter  ducts,  Fig.  5  shows 
several  constant  phase  lines  starting  out  of  the  lower  limit  of  the  search  region  and 
returning  to  the  lower-half  complex  qj  j/tj  plane.  Fig.  6  shows  the  presence  of  such  lines 
for  the  30  meter  and  the  40  meter  ducts.  These  constant  phase  lines  have  not  been 
observed  in  those  cases  of  lower  duct  heights. 
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Figure  5.  Upward  going  constant  phase  lines  in  the  qj  jAj  plane  for  the  10  meter  (top) 
and  20  meter  (bottom)  ducts. 
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Figure  6.  Upward  going  constant  phase  lines  in  the  q^/tj  plane  for  the  30  meter  (top) 


and  40  meter  (bottom)  ducts. 
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m.  ANALYSIS  AND  CONCLUSIONS 

A.      MODE  SEARCH  PARAMETERS 

The  implementation  of  a  mode  search  procedure  requires  several  parameters.  The 
search  region  has  to  be  defined  first.  In  M-Layer,  the  desired  maximum  range  attenuation 
rate  of  the  modes  which  support  wave  propagation  is  treated  as  an  input  parameter  called 
aloss'  S*ven  conventionally  in  dB  per  km.  This  parameter  determines  the  region  over 
which  the  modes  are  to  be  searched  as  explained  in  Chapter  I.  In  the  complex  q  plane, 
under  the  assumption  that  both  m^  and  p/k  of  interest  are  close  to  unity,  the  search  region 
is  bounded  from  above  by  the  FORTRAN  variable  qt      given  by 

q     .-*£*-.,      .  (12) 


The  location  at  which  to  start  the  mode  search,  qtest,  is  another  parameter 
determined  by  M-Layer.  Based  on  optical  considerations,  the  M-deficiency,  dmm|n  of 
Fig.  1,  which  is  the  amount  of  decrease  of  the  modified  refractivity  from  its  value  on  the 
sea  surface  to  its  minimum  value  in  the  air,  is  a  measure  of  the  capability  of  the  duct  to 
trap  electro-magnetic  waves.  As  explained  also  in  Chapter  I,  this  quantity  determines  the 
location  for  the  program  to  start  searching  for  modes.  The  validity  of  the  optical  argument 
and  the  usefulness  of  Eq.  (10)  have  been  proved  by  this  work,  as  pointed  out  in  Chapter 
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n. 

To  search  for  zeros  of  the  mode  function,  M-Layer  first  covers  the  search  region 
with  identical  mesh  squares.  It  then  follows  the  lines  of  the  type  Im{D(q)}  =  0  through 
each  mesh  square,  checking  for  indications  that  a  curve  Re{D(q)}  =  0  is  entering  the 
same  mesh  so  that  an  intersection  is  possible.  The  term  "mesh  size"  will  denote  the  size 
of  the  edges  of  the  mesh  squares.  It  is  the  size  of  the  step  taken  by  the  program  to 
advance  through  the  search  region.  It  also  determines  the  initial  resolution  of  the  location 
of  a  mode.  The  choice  of  the  mesh  size  should  strike  a  balance  between  the  desire  for  a 
speedy  completion  of  the  search  process  and  the  requirement  in  mode  locating  accuracy. 
As  reported  in  Ref.  2,  for  all  the  evaporation  ducts  considered,  the  choice  of  a  mesh  size 
equivalent  to  an  attenuation  rate  of  1/32  dB  per  km  appears  to  be  optimal,  that  is,  all 
modes  can  be  found  for  all  cases  investigated  in  Ref.  2  without  experiencing 
extraordinarily  long  computation  time.  Consider  this  mesh  size  as  the  default  and  call  it 
q^  in  the  complex  q  plane.  Under  the  same  assumption  for  deducing  Eq.  (12),  the  value 
of  q^/t^  at  9.6  GHz  equals  3.58x10  °.  Figure  7  shows  the  separation  between  two 
adjacent  constant  phase  lines  of  the  type  Im{D(q)}  =  0,  indicated  in  the  figure  as  solid 
lines,  along  the  top  edge  of  the  search  region  for  the  20  meter  duct.  The  minimum  of 
these  separations  in  terms  of  q^/tj  is  about  2.8x10  ,  which  corresponds  to  a  spacing 
of  a  little  less  than  eight  mesh  squares  apart.  The  minimum  separation  between  adjacent 
Im{D(q)}  =  0  and  Re{D(q)}  ■  0  lines  is  thus  about  four  meshes.  Constant  phase  line 
separations  for  all  evaporation  ducts  considered  are  included  in  Appendix  E.  Note  that  the 
minimum  separation  is  almost  constant  for  ducts  higher  than  20  meters.  It  increases 


22 


(U/l  lb)3^  J°  33U0J3JJIQ 


Figure  7.  Separation  of  Im{D(q)}  =  0  lines  along  the  top  edge  of  the  search  region 
(20  m  duct). 
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slightly  as  duct  height  is  decreased. 

As  the  program  follows  an  Im{D(q)}  =  0  line  into  the  search  region,  a  neighboring 
Re{D(q)}  =  0  line  moves  closer  and  enters  into  the  same  mesh  square.  This  causes  the 
program  to  invoke  the  root  finding  routine  to  determine  the  possible  locations  of  zeros 
of  the  mode  function  within  this  mesh  square.  If  the  mesh  size  is  too  large  and  the 
Re{D(q)}  =  0  line  enters  the  mesh  square  long  before  the  actual  intersection  takes  place, 
the  root  finding  routine  will  be  prematurely  invoked  many  times  and  the  probability  of 
producing  false  modes  is  increased.  This  explains  why  a  larger  mesh  size  sometimes  will 
increase  the  execution  time. 

Given  these  three  parameters  qt0  ,  qtesr  and  q^,  a  mode  search  strategy  which  does 
without  the  "contour  rectangles"  is  proposed.  It  should  improve  the  efficiency  of  M-Layer. 

B.      MODE  SEARCH  STRATEGY 

From  the  constant  phase  line  plots  of  Appendix  D,  a  strategy  to  search  for  the  mode 
can  be  drawn:  move  along  the  top  edge  starting  at  qtest>  Search  first  toward  either  the  left 
or  the  right,  then  reverse  course  to  search  along  the  other  direction.  After  the  search  along 
the  top  edge  is  completed,  search  the  lower  edge  from  one  end  to  the  other.  In  what 
follows,  the  implementation  of  this  strategy  will  be  discussed. 

1.       Track  Termination 

The  tracking  of  an  Im{D(q)}  =  0  line  naturally  ends  if  the  top  edge  or  the 
bottom  edge  of  the  search  region  is  reached.  On  the  other  hand,  the  search  region  is 
unbounded  to  the  right  and  left  of  qtest-  It  is  convenient  to  retain  the  feature  in  the 
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original  program  to  set  a  limit  on  the  number  of  steps  allowed  to  follow  a  constant  phase 
line.  From  Fig.  3,  this  limit  can  be  set  to  2.5  times  the  number  of  mesh  squares  between 
the  vertical  limits  of  the  search  region.  This  number  equals  2.5x32xajQSS. 

2.       Search  Termination 

M-Layer  has  to  determine  that  no  more  modes  within  the  search  region  is  to 
be  found  and  to  terminate  the  search  for  modes.  When  searching  along  the  top  edge  of 
the  search  region  for  the  constant  phase  lines  Im{D(q)}  =  0,  the  separation  between  the 
real  parts  of  the  mode  eigenvalues  is  of  interest  Figure  8  shows  the  separation  in  the  real 
part  of  neighboring  q-eigenvalues  scaled  by  tj,  plotted  against  their  locations  q^tj  along 
the  real  q^  j/tj  axis  for  the  case  of  the  20  m  duct.  Similar  plots  for  all  duct  heights  are 
grouped  in  Appendix  F.  Excluding  the  lowest  point  which  involves  the  mode  obtained  via 
searching  along  the  lower  edge,  the  separation  between  neighboring  modes  is  erratic  with 
an  upward  trend  away  from  the  center  of  the  figure,  which  is  close  to  the  search  starting 
position. 

The  increase  in  distance  between  two  modes  towards  the  end  of  the  range 
searched  makes  it  difficult  to  implement  an  adaptive  mechanism  to  terminate  the  search. 
On  the  other  hand,  Fig.  3  shows  that  almost  all  phase  lines  entering  the  search  region 
from  the  top  that  contain  a  mode  within  this  region  are  bunched  together.  Therefore,  the 
parameter  nst  set  equal  to  four  and  may  be  adjusted,  can  be  used  to  stop  the  search 
along  the  top  edge  when  four  consecutive  constant  phase  lines  tracked  turn  up  no  mode. 

The  search  along  the  real  q  axis  can  be  confined  to  within  the  end  points  of 
the  search  along  the  top  edge. 
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Figure  8.  Difference  between  consecutive  ReCq^j)  values  (20  m  duct). 
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3.       Track  Duplication  Avoidance 

When  the  program  searches  step  by  step  along  the  top  or  bottom  edges  for 
sign  changes  in  Im{D(q)}  to  start  tracking  the  constant  phase  line,  the  exit  of  a  constant 
phase  line  from  a  previous  track  will  be  encountered.  Entering  into  the  search  region  at 
such  a  location  will  simply  retrace  a  constant  phase  line  which  has  already  been  examined 
for  the  existence  of  a  mode.  Thus,  the  exit  locations  of  the  constant  phase  lines  must  be 
recorded  and  checked  whenever  a  sign  change  in  Im{D(q)}  is  found  before  deciding  to 
follow  the  constant  phase  line  into  the  search  region.  Along  the  top  edge,  a  record  of  four 
updated  most  recent  exit  locations  from  the  top  edge  should  be  adequate.  The  locations 
of  exits  from  the  bottom  should  be  kept  for  use  during  the  search  along  the  lower  edge 
of  the  search  region.  For  this  record,  a  dimension  of  256  should  be  adequate  for  the 
current  cases.  But  this  dimension  should  be  adjusted  if  other  types  of  ducts  are  studied. 
A  record  of  four  of  the  most  recent  exit  locations  out  of  the  lower  edge  should  also  be 
kept  and  checked  to  avoid  duplicate  tracking  of  the  constant  phase  lines. 
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APPENDIX  A:  MODE  LOCATIONS  IN  THE  COMPLEX  qil/tj  PLANE 

This  appendix  contains  figures  of  mode  locations  of  evaporation  ducts  of  2,  4,  6, 
10,  20,  30,  and  40  m  heights  plotted  in  the  complex  q^/tj  plane  for  easy  comparison. 
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Figure  A.l   q^/tj  mode  locations  (2  m  duct). 
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Figure  A.2  q^/tj  mode  locations  (4  m  duct). 
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Figure  A.3  q^jAj  mode  locations  (6  m  duct). 
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Figure  A.5   qji/tj  mode  locations  (10  m  duct). 
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Figure  A.6  <\\\h\  mode  locations  (12  m  duct). 
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Figure  A.7  qij/tj  mode  locations  (14  m  duct). 
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Figure  A.8  qij/tj  mode  locations  (16  m  duct). 
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Figure  A.9  qj^/tj  mode  locations  (18  m  duct). 
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Figure  A.10  q^Aj  mode  locations  (20  m  duct). 
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Figure  A.  12  qjj/tj  mode  locations  (24  m  duct). 
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Figure  A.  13  qij/tj  mode  locations  (26  m  duct). 
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Figure  A.14  qij/tj  mode  locations  (28  m  duct). 
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Figure  A.15  qjjAj  mode  locations  (30  m  duct). 
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Figure  A.16  q^jAj  mode  locations  (32  m  duct). 
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Figure  A.17  qji/tj  mode  locations  (34  m  duct). 
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Figure  A.18   qjj/tj  mode  locations  (36  m  duct). 
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Figure  A.19  q ^  j/t ^  mode  locations  (38  m  duct). 
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Figure  A.20  qnAj  mode  locations  (40  m  duct). 
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APPENDIX  B:  SUBROUTINE  FNDMOD  AND  FZEROX 

This  appendix  contains  the  listing  of  the  subroutines  FNDMOD  and  FZEROX, 
modified  from  the  NPS  version  of  the  M-Layer  FORTRAN  code  to  track  the  constant 
phase-lines  and  to  locate  the  modes. 
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Line#  Source  Line       Microsoft  FORTRAN  Optimizing  Compiler  Version  5.00 
1  subroutine  fndmod(aloss,dmmin,tl,dmesh,filem4umode,qeigen) 

2 

3  c  purpose: 

4  c       This  subroutine  sets  up  the  areas  in  the  complex  qll-plane 

5  c       to  search  for  the  zeroes  of  the  modal  function. 

6  c 

7  c  inputs... 

8  c  mxlayr  -  maximum  number  of  layers  allowed  in  refractivity 

9  c  profile 

10  c  nzlayr  -  actual  number  of  layers  in  refractivity  profile 

1 1  c  aloss  -  maximum  attenuation  rate  (in  db/km)  of  modes 

12  c  to  be  found 

13  c  dmmin  -  minimum  of  zim(j)-zim(l) 

14  c  tl       -koal23 

15  c  dmesh  -  initial  mesh  size  divisor 

16  c 

17  c  outputs... 

18  c  qeigen  -  complex  array  containing  the  zeros  of  the  modal 

19  c  function 

20  c  nrmode  -  actual  number  of  modes  found 

21  c 

22  c  subroutines  called... 

23  c  fzerox 

24  c  findfx 

25  c  common  block  areas... 

26  c  coml 

27  £********** 

28  c 

29  implicit  double  precision  (a-h,o-z) 

30  complex*  16  ctemp,tl,qeigen,zeros,fxl,fx2 

3 1  parameter(c201og=8.68588963806504d0,one=(  1  .d0,0.d0), 

32  $  Step0=1024.d0,step0h=step0/2.d0) 

33  character*3     filemjdine 

34  character*40  rfile 

35  c 

37  c       qeigen     -  complex  array  containing  all  the  zeros  of  the 

38  c  modal  function  found 

39  c       zeros       -  complex  array  containing  the  zeros  of  the  modal 

40  c  function  found  in  the  current  rectangular  region 

41  c  of  the  complex  qll-plane 

42  c 

44  c         use  include  file  for  parameters  of 

45  c  mxlayr  max  #  layers 

46  c  mxmode  max  #  modes 

47  c 

48  $include:  'mlaparm.inc' 

*****  Begin  listing  of:  mlaparm.inc 
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1  c 

2  c         include  file  to  define  the 

3  c  maximum  #  of  layers  (mxlayr) 

4  c  maximum  #  of  modes  (mxmode) 

5  c 

6  parameter  (mxlayr=35  ) 

7  parameter  (mxmode=127  ) 
*****  End  listing  of:  mlaparm.inc 

49 

50  dimension  qeigen(mxmode),zeros(2*mxmode+l) 

51  c***** 

52  common /com  1/waveno 

53  c 

55  c  rl      -  value  of  mode  search  on  the  real  part  of  ql  1  at  the 

56  c  left  edge  in  tmesh  units. 

57  c  r2     -  value  of  mode  search  on  the  real  part  of  ql  1  at  the 

58  c  right  edge  in  tmesh  units. 

59  c  bot  -  value  of  the  imaginary  part  of  ql  1  at  the  bottom 

60  c  edge  (this  is  set  to  0). 

61c       tO      -  value  of  the  imaginary  part  of  ql  1  at  the  top  edge 

62  c  in  tmesh  units. 

63  c       step  -  size  of  search  areas. 

65  c 

66  c 

67  c 

68  c  start  of  executable  statements 

69  c 

70  c* *********************************************************************** 

71c       set  up  search  areas  for  finding  modes  and  solve  for  modes, 

72  c       calculate  approximate  value  for  re(qll)  where  modes  start 

73  c       occurring 

75  c 


76 

nrmode=0 

77  c 

78 

recons=dreal(tl) 

79 

rstart=-2.0d-6*dmmin 

80 

qtest=rstart*recons 

81 

82 

ttopl=(2.0d-3/(waveno*c201og))*recons 

83 

ttop=aloss*ttopl 

84 

tmesh=ttop  1/dmesh 

85 

if  (tmesh  .gt.  0.1  dO)  tmesh=1.0d-l 

86 

if(tmesh  .gt.  1.0d-3)  then 

87 

tol=1.0d-4 

88 

else 

89 

tol=tmesh*0.1d0 

90 

end  if 
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91   c 

92  c***** 

93 

write(M002) 

94 

write(16,1002) 

95  c 

96 

write(M003)  tmesh,tol 

97 

write(16,1003)  tmesh.tol 

98  c 

99 

tO=dnint(ttop/tmesh)+ 1  .d0 

100 

bot=0.d0 

101   c** 

*** 

102  c 

103 

rl=dnint(qtest/tmesh)-stepOh 

104 

rl=rl 

105 

rr=rl+step0 

106 

call  findfx(rl,tO,fxl,xl,yl,tmesh) 

107 

il=int(rl) 

108 

in=int(step0)+int(rl) 

109 

nline=0 

110  50 

continue 

111 

do  100  nn=il,in 

112 

r2=rl+l.d0 

113 

call  findfx(r2,t0,fx2,x2,y2,tmesh) 

114 

if  (yl*y2  .It.  O.dO)  then 

115 

r=rl*tmesh 

116 

write  (*,5555)  r 

117  c 

write  (16,5555)  r 

118  5555       format  (/5x,'r=  '428.16) 

119 

go  to  400 

120 

endif 

121 

rl=r2 

122 

fxl=fx2 

123 

xl=x2 

124 

yl=y2 

125    100    continue 

126 

return 

127  c 

128  400     continue 

129 

nrold=nrmode 

130 

nline=nline+l 

131 

ncl00=int(nline/100) 

132 

nc  10=int((nline-nc  100*  100)/10) 

133 

nc  l=nline-nc  100*  100-nc  10*  10 

134 

ktine=char(48+ncl00)//char(48+ncl0)//char(48+ncl) 

135 

rfile='e:Nans5V7/filem//kline//'.dat' 

136  c****** 

137 

open  (unit=25,file=rfile,status='unknown') 

138 

139 

call  fzerox(t0j-l /1 1  ,fx  1  ,x  1  ,y  1  ,fx2,x2,y2,zeros,nrold,m 

140 

$               tmesh) 
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141  close  (25) 

142  c****** 

143  c 

144  nrmode=nrnew 

145  newO=nrmode-nrold 

146  c 

147  r* ************************************************** ***************** 

148  c       mode  search  completed 

149  c       order  zeros  found  by  order  of  increasing  real  part. 

150  q**** *********************************************************  ******* 

151  c 

152  if(nrmode  .gt.  1)  then 

153  jkend=nrmode- 1 
154 

155  do420jk=ljkend 

156  nend=nrmode-jk 
157 

158  do410ja=l,nend 

159  nrj=nrmode-ja 

160  nrjl=nrj+l 

161  if(dimag(zeros(nrjl)).lt. 

162  $  dimag(zeros(nrj)))  then 

1 63  ctemp=zeros(nrj  1 ) 

164  zeros(nrjl)=zeros(nrj) 

165  zeros(nrj)=ctemp 

166  end  if 

167  410         continue 

168  420       continue 

169  end  if 

170  c 

171  ^* ********************************************************************* 

172  c    the  possibility  exists  that  duplicate  (within  the  tolerance  'tol') 

173  c    zeros  of  the  modal  equation  will  be  found,  eliminate  these 

174  c    duplicate  zeros. 

175  c********************************************************************** 

176  c 

177  jkflag=0 

178  jkend=nrmode-l 
179 

180  do240jk=ljkend 

181  jkl=jk+l 

182  ctemp=zeros(jkl) 

183  chksq=cdabs(zeros(jk-jkflag)-ctemp) 

184  if(chksq  .It.  tol)  then 

185  jkflag=jkflag+l 

186  go  to  240 

187  end  if 

188  zeros(jkl-jkfIag)=ctemp 

189  240    continue 

190  nrmode=nrmode-jkflag 
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191  c 

192  c* ************************** **************** 

193  c    Store  the  zeros  as  the  eigenvalues. 

j 94   c********** ********* ************* *********** 

195  c 

196  nrmode=minO(nrmode,mxmode) 
197 

198  do  430  jk=  1  jirmode 

199  qeigen(jk)=zeros(jk) 

200  430    continue 

201  c 

202  il=int(rll)+l 

203  if  (il  .It.  in)  then 

204  rl=rll+hd0 

205  caU  findfx(rl,t0,fxl,xl,yl,tmesh) 

206  go  to  50 

207  endif 

208  return 

209  c 

210  c 

211  c  format  statements 

212  c 

213  1000  format(/5x,'searching  for  zeroes  in  this  areas  are', 

214  $    '  defined  by:75x,'istart=  'jl0/5x,'itop=  \il0, 

215  $    5x,5  ibot=  ',ill/5x/ileft=  \il0,5x,'iright=  ',il0) 
216 

217   1001   format(5xj4,'  new  zeroes  found  in  this  area.'/) 
218 

219  1002  format(//5x,'*******        start  search  for  modal  eigenvalues', 

220  $         '        *******') 

221 

222   1003   format(/5x,'tmesh=  ',dl5.5,5x,'tol=  \dl5.5/) 

223 

224  end 

225  c 

226  c 

227  c* ***************************************** **************************** 

228  c* ********************************************************************* 

229  c 

230  subroutine  fzerox(t0,rl,rll,fxl,xl,yl,fx2,x2,y2,zeros,nrold, 

231  $  nrnew,tmsh0) 

232  c 

233  c***** 

234  c  fzerox  is  a  routine  for  finding  the  zeroes  of  a  complex  function,  f, 

235  c  which  lie  within  a  specified  rectangular  region  of  the 

236  c  complex  qll  plane,  assuming  that  the  function  has  only 

237  c  simple  zeroes  over  this  rectangle. 

238  c 

239  c  parameters  specifying  the  search  rectangle: 

240  c  tmesh  -  set  equal  to  about  half  the  average  spacing  between 
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241  c  zeroes  within  the  rectangle.  A  smaller  value  may  be  used 

242  c  as  a  safety  measure,  but  too  small  a  value  will  result 

243  c  in  excessively  long  run  time. 

244  c       zeros  -  output  list  of  (complex)  values  of  ql  1  at  which 

245  c  zeroes  are  found. 

246  c  nrnew.nrold-  the  number  of  zeroes  found 

247  c 

248  c  subroutines  called- 

249  c  findfx 

250  c  roots 

251  c***** 

252  implicit  double  precision  (a-h,o-z) 

253  complex*16  fl0,f01,fll,rxl,fx2,fx00ixl0ix01,fxll, 

254  +  one,sol,zeros 

255  parameter(one=(  1  .d0,0.d0)) 

256  $include:  'mlaparm.inc' 

*****  Begin  listing  of:  mlaparm.inc 

1  c 

2  c         include  file  to  define  the 

3  c  maximum  #  of  layers  (mxlayr) 

4  c  maximum  #  of  modes  (mxmode) 

5  c 

6  parameter  (mxlayr=35  ) 

7  parameter  (mxmode=127  ) 
*****  End  listing  of:  mlaparm.inc 

257  dimension  sol(3),theta(2),zeros(2*mxmode+l),rline(2,1024) 

258  c 

259  common/tmccom/tmesh 

260  c***** 

261  c 

262  npoint=0 

263  nrnew=nrold 

264  tmesh=tmshO 

265  bot=0.d0 

266  t=t0-l.d0 

267  rll=rl 

268  r=rll 

269  c***** 

270  c 

271  fx01=fxl 

272  fx01r=xl 

273  fx01i=yl 

274  fxll=fx2 

275  fxllr=x2 

276  fxlli=y2 

277  go  to  82 

278  c 

279  c***** 

280  c  enter  mesh  square  from  left  side  or  exit  rectangle  at  right  edge. 

281  c 
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282  20 

r=r+l.d0 

283   21 

fx01=fxll 

284 

fx01r=fxllr 

285 

fx01i=fxlli 

286 

fx00=fxl0 

287 

fx00r=fxl0r 

288 

fx00i=fxl0i 

289  22 

continue 

290 

call  findfx(r+ 1  .d0,t+ 1  .d0,fx  1 1  ,fx  1  lrix  1 1  i,tmesh) 

291 

call  findfx(r+ 1  .d0,Ufxl0 Jx  10r,fxl0i,tmesh) 

292  q******* 

293  c 

Determine  the  edge  of  exit  of  im(f)=0  from  current  mesh. 

294 

edgeit=fx01i*fxlli 

295 

edgeib=fx00i*fxl0i 

296 

if  (edgeib  .gt.  O.dO)  then 

297   c 

lm(f)=0  goes  through  the  01  to  10  line. 

298 

if  (edgeit  .gt.  O.dO)  then 

299  c 

lm(f)=0  goes  through  the  10  to  11  edge  (edge  1). 

300 

lout=l 

301 

else 

302  c 

lm(f)=0  goes  through  the  01  to  11  edge  (edge  2) 

303 

lout=2 

304 

end  if 

305 

else 

306  c 

lm(f)=0  goes  through  the  00  to  10  edge  (edge  4) 

307 

lout=4 

308 

if  (edgeit  .lL  O.dO)  then 

309  c 

lm(f)=0  also  runs  through  01  to  11  and  10  to  11  edges. 

310  c 

Store  crossing  location  and  in/out  information. 

311 

knot34=knot34+l 

312  c 

loc34r(knot34)=lr 

313  c 

loc34i(knot34)=li 

314 

end  if 

315 

end  if 

316  c******* 

317 

go  to  85 

318  c***** 

319  c  enter  mesh  square  from  bottom  side  or  exit  rectangle  at  top  edge. 

320  40 

t=t+l.d0 

321  41 

fx00=fx01 

322 

fx00r=fx01r 

323 

fx00i=fx01i 

324 

fxlO=fxll 

325 

fxl0r=fxllr 

326 

fxl0i=fxlli 

327  42 

continue 

328 

callfmdfx(r,t+l.d0/x01ix01r,fx01i,tmesh) 

329 

call  findfx(r+ 1  .d0,t+ 1  .d0,fx  1 1  ,fx  1 1  rix  1 1  i.tmesh) 

330  c******* 

331   c 

Determine  the  edge  of  exit  of  im(f)=0  from  current  mesh. 
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332 

edgeil=fx00i*fx01i 

333 

edgeir=fxlOi*fxlli 

334 

if  (edgeir  .gt.  O.dO)  then 

335   c 

lm(f)=0  goes  through  the  00  to  1 1  line. 

336 

if  (edgeil  .gt.  O.dO)  then 

337  c 

lm(f)=0  goes  through  the  01  to  11  edge  (edge  2) 

338 

lout=2 

339 

else 

340  c 

lm(f)=0  goes  through  the  00  to  01  edge  (edge  3). 

341 

lout=3 

342 

end  if 

343 

else 

344   c 

lm(f)=0  goes  through  the  10  to  11  edge  (edge  1) 

345 

lout=l 

346 

if  (edgeil  .It.  O.dO)  then 

347   c 

lm(f)=0  also  runs  through  00  to  01  and  01  to  11  edges. 

348  c 

Store  crossing  location  and  in/out  information. 

349 

knot41=knot41+l 

350  c 

loc41r(knot41)=lr 

351   c 

Ioc41i(knot41)=li 

352 

end  if 

353 

end  if 

354  c*< 

<***** 

355 

go  to  85 

356  c*< 

>*** 

357   c  < 

inter  mesh  square  from  right  side  or  exit  rectangle  at  left  edge. 

358 

359  60 

r=r-l.d0 

360  61 

fxll=fx01 

361 

fxllr=fx01r 

362 

fxlli=fx01i 

363 

fxl0=fx00 

364 

fxl0r=fx00r 

365 

fxl0i=fx00i 

366  65 

continue 

367 

call  findfx(r,t+ 1  .dOixO  1  ,f xO  lr,fx0 1  Umesh) 

368 

call  findfx(r.t,fx00,fx00rix00i,tmesh) 

369  c*< 

>***** 

370  c 

Determine  the  edge  of  exit  of  im(f)=0  from  current  mesh. 

371 

edgeit=fx01i*fxlli 

372 

edgeib=fx00i*fxl0i 

373 

if  (edgeit  .gt.  O.dO)  then 

374  c 

lm(f)=0  goes  through  the  01  to  10  line. 

375 

if  (edgeib  .gt.  O.dO)  then 

376  c 

lm(f)=0  goes  through  the  00  to  01  edge  (edge  3). 

377 

lout=3 

378 

else 

379  c 

Im(f)=0  goes  through  the  00  to  10  edge  (edge  4) 

380 

lout=4 

381 

end  if 
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382 

else 

383  c 

lm(f)=0  goes  through  the  01  to  11  edge  (edge  2) 

384 

lout=2 

385 

if  (edgeib  .It.  O.dO)  then 

386  c 

lm(f)=0  also  runs  through  00  to  10  and  00  to  01  edges. 

387  c 

Store  crossing  location  and  in/out  information. 

388 

knot  12=knot  12+1 

389  c 

locl2r(knotl2)=lr 

390  c 

locl2i(knotl2)=li 

391 

end  if 

392 

end  if 

393  c** 

***** 

394 

go  to  85 

395  c***** 

396  c  enter  mesh  square  top  side  or  exit  rectangle  at  bottom  edge. 

397  c 

398  80 

t=t-l.d0 

399  81 

fx01=fx00 

400 

fx01r=fx00r 

401 

fx01i=fx00i 

402 

fxll=fxl0 

403 

fxllr=fxl0r 

404 

fxlli=fxl0i 

405   82 

continue 

406 

call  findfx(r,t,fx00,fx00r,fx00i,tmesh) 

407 

caUfindfx(r+l.d0,t,fxl0Jxl0r,fxl0i,tmesh) 

408  c 

409  c******* 

410  c 

Determine  the  edge  of  exit  of  im(f)=0  from  current  mesh. 

411   c 

412   83 

continue 

413 

edgeil=fx00i*fx01i 

414 

edgeir=fxl0i*fxlli 

415 

if  (edgeil  .gt.  O.dO)  then 

416  c 

lm(f)=0  goes  through  the  00  to  1 1  line. 

417 

if  (edgeir  .gt.  O.dO)  then 

418  c 

lm(f)=0  goes  through  the  00  to  10  edge  (edge  4) 

419 

lout=4 

420 

else 

421   c 

lm(f)=0  goes  through  the  10  to  11  edge  (edge  1). 

422 

lout=l 

423 

end  if 

424 

else 

425  c 

lm(f)=0  goes  through  the  00  to  01  edge  (edge  3) 

426 

lout=3 

427 

if  (edgeir  .It.  O.dO)  then 

428  c 

lm(f)=0  also  runs  through  00  to  10  and  10  to  1 1  edges. 

429  c 

Store  crossing  location  and  in/out  information. 

430 

knot23=knot23+l 

431   c 

loc23r(knot23)=lr 
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432  c  loc23i(knot23)=li 

433  end  if 

434  end  if 

435  c************************************************************************ 

436  c     Test  for  there  being  at  least  one  re(f)=0  line  entering  and 

437  c     leaving  the  mesh  square. 

438  c************************************************************************ 

439  c 

440  85      continue 

441  npoint=npoint+l 

442  rrr=r*tmesh 

443  ttt=t*tmesh 

444  rline(ljipoint)=rrr 

445  rline(2,npoint)=ttt 

446  if  (npoint  .eq.  1024)  lout=5 

447  c 

448  if((t  .It.  bot)  .or.  (t  .gt.  t0))  go  to  100 
449 

450  if  ((fx00r*fxl0r  .gt.  O.dO)  .and.  (fx01r*fxl lr  .gt.  O.dO) 

451  +   .and.  (fx00r*fx01r  .gt.  O.dO))  go  to  (20,40,60,80,100) 

452  +  lout 

453  q************************************************************************* 

454  c      Computate  the  values  of  the  modal  function  at  the  corners  of 

455  c      a  mesh  square  to  determine  its  Taylor  series  to  the  3rd  order 

456  c      for  estimating  its  root  locations. 

457  c****************** ************************** ************************** 

458  c 

459  c       f00=one 

460  fl0=cdexp(fxl0-fx00)-one 

461  f01=cdexp(fx01-fx00)-one 

462  f  1  l=cdexp(fx  1  l-fx00)-one 

463  c 

464  c***********  estimate  locations  of  zeroes  by  radicals***************** 

465  c 

466  caUroots(flOJOUlUoljirsol) 

467  c 

468  do  90  n=l,nrsol 

469  ureal  =  dreal(sol(n)) 

470  uimag  =  dimag(sol(n)) 

471  if  (ureal  .It.  O.dO  .or.  ureal  .gt.  l.OdO)  go  to  90 

472  if  (uimag  .It.  O.dO  .or.  uimag  .gt.  l.OdO)  go  to  90 

473  92        theta(l)=(r+ureal)*tmesh 

474  theta(2)=(t+uimag)*tmesh 

475  nrnew  =  nrnew+1 

476  zeros(nrnew)=dcmplx(theta(l),theta(2)) 

477  90     continue 

478  c 

479  c* ****************************************** 

480  c      continue  following  the  phase  line 

481  c* ****************************************** 
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482 

c 

483 

95 

continue 

484 

go  to  (20,40,60,80,100)  lout 

485 

c** 

486 

100 

continue 

487 

Q****** 

488 

write  (25,8888)  ((rline(ij),  i=l, 

2),j=l, 

,  npoint) 

489 

8888  f 

490 
491 
492 

c 

return 
end 
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APPENDIX  C:  CONSTANT  PHASE  LINE  TRACKING  FLOW  CHARTS 


START 


nrmode  -  0 
recons  =  dreal(t1) 
qtest=-2x10A(-6)  x  dmmin  x  recons 
ttop1=((2x10A(-3)  /  k  x  20log(e))  x 

recons 
top-  alossxttopl 
tmesh  -  ttopl  /  dmesh 


tol  -  tmesh  x  0.1 


mesh  =  0.1 


tol  -  0.0001 
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tO=dnint(ttop/tmesh)+1  .dO 

bot  -  O.dO 

r1  =dnint(qtest/tmesh)-stepOh 

rl-r1 

rr-r1+1024.d0 


CALL    FINDFX 


i1  -  int(M)  ,in  -  1024+int(r1),  nline=0 


DO  100 
nn=M,in 
r2=M+1.dO 


CALL    FINDFX 


r=M  x  tmesh 

/  WRITE  ( r )  I 
(400) 
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MOC^          w 

KZ 

r           i 

nrold  -  nrmode 

nline  -  nline  +  1 

nc100  =  int(nline/100) 

nd  0   -  int((nline-nd  00*1 00)/1 0) 

nd    -  nline  -  nd  00*1 00  -  nd  0*1 0 
kline=char(48+nd  00)//char(48+nc1 0 

)//char(48+nd) 
rfile-,e:\mlxhe\ans\r7/filem//kline//.dat' 

4r 

Open(unit=25,  file=rfile,  status= 
'unknown') 

* 

CALL      FZEROX 

+ 

CLOSE  ( 25 ) 

nrmode  -  nmew 
newO  -  nrmode  -  nrold 
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WRITE 
(newO) 


K 


CONTINUE 


DO  420 
jk  ■  1  ,  jkend 
nend  -  nrmode  -  jk 


ja  -  1  ,  nend 
nrj  -  nrmode -ja 
nrjl  -  nrj  +  1 


ctemp  -  zeros  (nrjl) 
zeros  (niJ1)  -  zeros  (nrj) 
zeros  (nrj)  -  ctemp 


0 

nrm<xto-1 
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/ 


DO  240 

jk-1,jkend 


CONTINUE 


v        jki-|k+i 
\       ctemp- zeros  Qk1) 


chksq  - 
cdabs(zeros 
flk-jkflag)- 
ctemp) 


Jkflag-jkflag+1 


CONTINUE 


zeros()k1-jkflao)  -  ctemp 


nrmode=nrmode  - 

jkflag 
nrmode-minO(nrmode, 
mxmode) 

DO  430 

jk  =  1  ,  nrmod 


qeigen(jk)  -  zerosflk) 

i 


i1  -int(r11)  +  1 


CALL  FINDFX 
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APPENDIX  D:  CONSTANT  PHASE  LINES  IN  THE  COMPLEX  qn/t1  PLANE 

This  appendix  contains  plots  of  the  constant  phase  lines  Im{D(q)}  =  0  and 
Re{D(q)}  =  0  initiating  from  the  top  edge  of  the  search  region  in  the  complex  q^/tj 
plane  for  the  evaporation  ducts  of  2,  4,  6,  8,  10,  20,  30,  and  40  m  heights. 
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Figure  D.l   Constant  phase  lines  in  the  q^j/tj  plane  (2  m  duct). 
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Figure  D.2   Constant  phase  lines  in  the  qjj/tj  plane  (4  m  duct). 
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Figure  D.3   Constant  phase  lines  in  the  qjjAj  plane  (6  m  duct). 
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Figure  D.4  Constant  phase  lines  in  the  qj  j/tj  plane  (8  m  duct) 
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Figure  D.5  Constant  phase  lines  in  the  q^/tj  plane  (10  m  duct) 
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Figure  D.6  Constant  phase  lines  in  the  qj  j/t  j  plane  (20  m  duct) 
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Figure  D.7   Constant  phase  lines  in  the  qjj/tj  plane  (30  m  duct) 
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Figure  D.8   Constant  phase  lines  in  the  qj  j/tj  plane  (40  m  duct). 
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APPENDIX  E:  SEPARATION  OF  Im{D(q)}=0  LINES 

This  appendix  contains  figures  of  the  separation  of  Im{D(q)}=0  lines  along  the  top 
edge  of  the  search  region.  The  2,  4,  6,  8,  10,  20,  30  and  40  meter  evaporation  ducts  are 
included. 
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Figure  E.1   Separation  of  Im{D(q)}=0  lines  along  the  top  edge  of  the  search  region 
(2  m  duct). 
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Figure  E.2  Separation  of  Im{D(q)}=0  lines  along  the  top  edge  of  the  search  region 
(4  m  duct). 
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Figure  E.3  Separation  of  Im{D(q)}=0  lines  along  the  top  edge  of  the  search  region 
(6  m  duct). 


78 


Figure  E.4  Separation  of  Im{D(q)}=0  lines  along  the  top  edge  of  the  search  region 
(8  m  duct). 
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Figure  E.5   Separation  of  Im{D(q)}=0  lines  along  the  top  edge  of  the  search  region 
(10  m  duct). 
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Figure  E.6  Separation  of  Im{D(q)}=0  lines  along  the  top  edge  of  the  search  region 
(20  m  duct). 
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Figure  E.7   Separation  of  Im{D(q)}=0  lines  along  the  top  edge  of  the  search  region 
(30  m  duct). 
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Figure  E.8   Separation  of  Im{D(q)}=0  lines  along  the  top  edge  of  the  search  region 
(40  m  duct). 
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APPENDIX  F:  DIFFERENCE  BETWEEN  CONSECimVE  ReCq^tj)  VALUES 

Separation  in  real  part  of  neighboring  q-eigenvalues  scaled  by  tj  is  plotted  in  this 
ippendix  against  the  eigenvalue  locations  along  the  real  q^  j/tj  axis. 
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Figure  F.l   Difference  between  consecutive  ReCq^tj)  values  (2  m  duct). 
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Figure  F.2  Difference  between  consecutive  Re(q^ti)  values  (4  m  duct). 
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Figure  F.3  Difference  between  consecutive  ReCq^j)  values  (6  m  duct). 
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Figure  F.4  Difference  between  consecutive  Re(q^tt)  values  (8  m  duct). 
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Figure  F.5  Difference  between  consecutive  ReCq^j)  values  (10  m  duct). 
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Figure  F.6  Difference  between  consecutive  ReCq^tj)  values  (20  m  duct). 
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Figure  F.7  Difference  between  consecutive  ReCq^tj)  values  (30  m  duct). 


91 


Figure  F.8   Difference  between  consecutive  ReCq^tj)  values  (40  m  duct). 
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